

# - Class Fixed Effects only. This is the main factor I should be
# concerned with.
# - Including SES variables and other variables. I should have up to
# three models and compare then. One with individual charactetics,
# other with indiv. + SES and other controlling for class FE.

ind_controls <- paste("male", "factor(religion_simple)", 
                      "age_grade_distortion", "score_self_esteem", 
                      "score_parents_support", "score_study", sep = "+")

ses_controls <- paste("scores_poverty", "score_neighborhood_quality", 
                      "factor(occ_father_simple_no_job_no_na)",
                      "factor(occ_mother_simple_no_job_no_na)", sep = "+")

class_fe <- paste("factor(class_id)")


general_reg_separe_race_het <- function(dep_var, grade_var, hetero_var) {
  # dep_var = "ssi_same_negro_na_2_std"
  # grade_var = "grades_score"
  # hetero_var = "EM"
  
  basic_reg <- paste0(dep_var, "~", paste0(grade_var))
  
  Form_basic <- formula(basic_reg)
  
  # Individual control
  Form_ind <- formula(paste(basic_reg, ind_controls, sep = "+"))
  
  # SES control
  Form_ses <- formula(paste(basic_reg, ses_controls, sep = "+"))
  
  # Individual + SES control
  Form_ind_ses <- formula(paste(basic_reg, ind_controls, 
                                ses_controls, sep = "+"))
  
  # Individual controls + Classroom FE
  Form_ind_class_fe <- formula(paste(basic_reg, ind_controls, 
                                     class_fe, sep = "+"))
  
  # SES controls + Classroom FE
  Form_ses_class_fe <- formula(paste(basic_reg, ses_controls, 
                                     class_fe, sep = "+"))
  
  # Ind. + SES controls + Classroom FE
  Form_ind_ses_class_fe <- formula(paste(basic_reg, ind_controls,
                                         ses_controls, class_fe, 
                                         sep = "+"))
  
  list_specifications <- c(Form_basic, Form_ind, Form_ses, Form_ind_ses,
                           Form_ind_class_fe, Form_ses_class_fe, 
                           Form_ind_ses_class_fe) 
  
  results_final_0 <- list()
  data_0 <- data[which(data[,hetero_var]==0),]
    
    print(1)
    results <- lapply(list_specifications, function(x) (lm(x, data = data_0,
                                                           subset = (race == 1))))
    results_r <- lapply(results, robust_se, "HC1")
    results_c <- lapply(results, cluster_se, "class_id")
    #results_c <- NA
    
    results_final_0[[1]] <- list(results, results_r, results_c)  
  
  print(4)
  results <- lapply(list_specifications, function(x) (lm(x, data = data_0, 
                                                         subset = (negro_na == 1))))
  results_r <- lapply(results, robust_se, "HC1")
  results_c <- lapply(results, cluster_se, "class_id")
  #results_c <- NA
  
  
  results_final_0[[4]] <- list(results, results_r, results_c)  
  
  results_final_1 <- list()
  data_1 <- data[which(data[,hetero_var]==1),]
    
    print(1)
    results <- lapply(list_specifications, function(x) (lm(x, data = data_1, 
                                                           subset = (race == 1))))
    results_r <- lapply(results, robust_se, "HC1")
    results_c <- lapply(results, cluster_se, "class_id")
   # results_c <- NA
    
    results_final_1[[1]] <- list(results, results_r, results_c)  
  
  print(4)
  results <- lapply(list_specifications, function(x) (lm(x, data = data_1,
                                                         subset = (negro_na == 1))))
  results_r <- lapply(results, robust_se, "HC1")
  results_c <- lapply(results, cluster_se, "class_id")
  #results_c <- NA
  
  
  results_final_1[[4]] <- list(results, results_r, results_c)  
  return(list(results_final_0, results_final_1))
}

variables_to_omit <- c("(Constant)|(class_)|(occ_)|(parents)|(self)|(study)|
                       |(poverty)|(neighborhood)|(religion)|(age)|(male)")

# High School and Elementary School
sep_same_linear_em <- general_reg_separe_race_het("ssi_same_negro_na_2_std",
                                                  "grades_score", 
                                                  hetero_var = "high_school")

sep_other_linear_em <- 
  general_reg_separe_race_het("ssi_other_race_negro_na_2_std",
                              "grades_score",
                              hetero_var = "high_school")

sep_all_linear_em <- general_reg_separe_race_het("ssi_all_race_2_std",
                                                 "grades_score",
                                                 hetero_var = "high_school")

list_regs_c_0 <- list(sep_all_linear_em[[1]][[1]][[3]][[7]], 
                      sep_all_linear_em[[1]][[4]][[3]][[7]],
                      sep_same_linear_em[[1]][[1]][[3]][[7]], 
                      sep_same_linear_em[[1]][[4]][[3]][[7]],
                      sep_other_linear_em[[1]][[1]][[3]][[7]],   
                      sep_other_linear_em[[1]][[4]][[3]][[7]])

list_regs_c_1 <- list(sep_all_linear_em[[2]][[1]][[3]][[7]], 
                      sep_all_linear_em[[2]][[4]][[3]][[7]],
                      sep_same_linear_em[[2]][[1]][[3]][[7]],
                      sep_same_linear_em[[2]][[4]][[3]][[7]],
                      sep_other_linear_em[[2]][[1]][[3]][[7]],   
                      sep_other_linear_em[[2]][[4]][[3]][[7]])

list_regs_0 <- list(sep_all_linear_em[[1]][[1]][[1]][[7]], 
                    sep_all_linear_em[[1]][[4]][[1]][[7]],
                    sep_same_linear_em[[1]][[1]][[1]][[7]], 
                    sep_same_linear_em[[1]][[4]][[1]][[7]],
                    sep_other_linear_em[[1]][[1]][[1]][[7]],   
                    sep_other_linear_em[[1]][[4]][[1]][[7]])

list_regs_1 <- list(sep_all_linear_em[[2]][[1]][[1]][[7]], 
                    sep_all_linear_em[[2]][[4]][[1]][[7]],
                    sep_same_linear_em[[2]][[1]][[1]][[7]], 
                    sep_same_linear_em[[2]][[4]][[1]][[7]], 
                    sep_other_linear_em[[2]][[1]][[1]][[7]],   
                    sep_other_linear_em[[2]][[4]][[1]][[7]])

n_obs_0 <- sapply(list_regs_0, get_obs)
n_obs_1 <- sapply(list_regs_1, get_obs)

stargazer(list_regs_c_0,
          column.separate = c(2,2,2),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White",
                             "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_0)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_em_0.tex")

stargazer(list_regs_c_1,
          column.separate = c(2,2,2),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c( "All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White", "Nonwhite",
                             "White", "Nonwhite"),
                           c("Observations", n_obs_1)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_em_1.tex")

  # Gender differences
sep_same_linear_gender <- 
  general_reg_separe_race_het("ssi_same_negro_na_2_std", 
                              "grades_score", 
                              hetero_var = "male")
sep_other_linear_gender <- 
  general_reg_separe_race_het("ssi_other_race_negro_na_2_std",
                              "grades_score", 
                              hetero_var = "male")
sep_all_linear_gender <- 
  general_reg_separe_race_het("ssi_all_race_2_std", 
                              "grades_score",
                              hetero_var = "male")

list_regs_c_0 <- list(sep_all_linear_gender[[1]][[1]][[3]][[7]], 
                      sep_all_linear_gender[[1]][[4]][[3]][[7]],
                      sep_same_linear_gender[[1]][[1]][[3]][[7]],
                      sep_same_linear_gender[[1]][[4]][[3]][[7]], 
                      sep_other_linear_gender[[1]][[1]][[3]][[7]],  
                      sep_other_linear_gender[[1]][[4]][[3]][[7]])
                    
list_regs_c_1 <- list(sep_all_linear_gender[[2]][[1]][[3]][[7]],
                      sep_all_linear_gender[[2]][[4]][[3]][[7]],
                      sep_same_linear_gender[[2]][[1]][[3]][[7]], 
                      sep_same_linear_gender[[2]][[4]][[3]][[7]], 
                      sep_other_linear_gender[[2]][[1]][[3]][[7]],   
                      sep_other_linear_gender[[2]][[4]][[3]][[7]])

list_regs_0 <- list(sep_all_linear_gender[[1]][[1]][[1]][[7]],
                    sep_all_linear_gender[[1]][[4]][[1]][[7]],
                    sep_same_linear_gender[[1]][[1]][[1]][[7]], 
                    sep_same_linear_gender[[1]][[4]][[1]][[7]], 
                    sep_other_linear_gender[[1]][[1]][[1]][[7]], 
                    sep_other_linear_gender[[1]][[4]][[1]][[7]])

list_regs_1 <- list(sep_all_linear_gender[[2]][[1]][[1]][[7]],
                    sep_all_linear_gender[[2]][[4]][[1]][[7]],
                    sep_same_linear_gender[[2]][[1]][[1]][[7]], 
                    sep_same_linear_gender[[2]][[4]][[1]][[7]], 
                    sep_other_linear_gender[[2]][[1]][[1]][[7]],  
                    sep_other_linear_gender[[2]][[4]][[1]][[7]])

n_obs_0 <- sapply(list_regs_0, get_obs)
n_obs_1 <- sapply(list_regs_1, get_obs)

stargazer(list_regs_c_0,
          column.separate = c(2,2,2),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White",
                             "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_0)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_gender_0.tex")

stargazer(list_regs_c_1,
          column.separate = c(2,2,2),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White", "Nonwhite",
                             "White", "Nonwhite"),
                           c("Observations", n_obs_1)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_gender_1.tex")

  # Share of nonwhite individuals:
data$sh_nw_class_dummy <- as.numeric(data$sh_nw_class >=
                                       median(data$sh_nw_class))

sep_same_linear_sh_nw_class <- 
  general_reg_separe_race_het("ssi_same_negro_na_2_std",
                              "grades_score", 
                              hetero_var = "sh_nw_class_dummy")
sep_other_linear_sh_nw_class <- 
  general_reg_separe_race_het("ssi_other_race_negro_na_2_std", 
                              "grades_score",
                              hetero_var = "sh_nw_class_dummy")
sep_all_linear_sh_nw_class <- 
  general_reg_separe_race_het("ssi_all_race_2_std", 
                              "grades_score",
                              hetero_var = "sh_nw_class_dummy")

list_regs_c_0 <- list(sep_all_linear_sh_nw_class[[1]][[1]][[3]][[7]], 
                      sep_all_linear_sh_nw_class[[1]][[4]][[3]][[7]],
                      sep_same_linear_sh_nw_class[[1]][[1]][[3]][[7]], 
                      sep_same_linear_sh_nw_class[[1]][[4]][[3]][[7]], 
                    sep_other_linear_sh_nw_class[[1]][[1]][[3]][[7]],  
                    sep_other_linear_sh_nw_class[[1]][[4]][[3]][[7]])

list_regs_c_1 <- list(sep_all_linear_sh_nw_class[[2]][[1]][[3]][[7]], 
                      sep_all_linear_sh_nw_class[[2]][[4]][[3]][[7]],
                      sep_same_linear_sh_nw_class[[2]][[1]][[3]][[7]],
                      sep_same_linear_sh_nw_class[[2]][[4]][[3]][[7]], 
                    sep_other_linear_sh_nw_class[[2]][[1]][[3]][[7]], 
                    sep_other_linear_sh_nw_class[[2]][[4]][[3]][[7]])

list_regs_0 <- list(sep_all_linear_sh_nw_class[[1]][[1]][[1]][[7]],
                    sep_all_linear_sh_nw_class[[1]][[4]][[1]][[7]],
                    sep_same_linear_sh_nw_class[[1]][[1]][[1]][[7]], 
                    sep_same_linear_sh_nw_class[[1]][[4]][[1]][[7]], 
                    sep_other_linear_sh_nw_class[[1]][[1]][[1]][[7]],   
                    sep_other_linear_sh_nw_class[[1]][[4]][[1]][[7]])

list_regs_1 <- list(sep_all_linear_sh_nw_class[[2]][[1]][[1]][[7]],
                    sep_all_linear_sh_nw_class[[2]][[4]][[1]][[7]],
                    sep_same_linear_sh_nw_class[[2]][[1]][[1]][[7]], 
                    sep_same_linear_sh_nw_class[[2]][[4]][[1]][[7]], 
                    sep_other_linear_sh_nw_class[[2]][[1]][[1]][[7]],  
                    sep_other_linear_sh_nw_class[[2]][[4]][[1]][[7]])

n_obs_0 <- sapply(list_regs_0, get_obs)
n_obs_1 <- sapply(list_regs_1, get_obs)

stargazer(list_regs_c_0,
          column.separate = c(4,4,4),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White", 
                             "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_0)),
          font.size = "footnotesize", type = "text")

stargazer(list_regs_c_1,
          column.separate = c(4,4,4),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite",
                             "White", "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_1)),
          font.size = "footnotesize", type = "text")

stargazer(list_regs_c_0,
          column.separate = c(4,4,4),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White", 
                             "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_0)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_sh_nw_class_0.tex")

stargazer(list_regs_c_1,
          column.separate = c(4,4,4),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White",
                             "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_1)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_sh_nw_class_1.tex")

# High and low grades of nonwhite individuals:
data$grades_score_dummy <- as.numeric(data$grades_score > 
                                        median(data$grades_score, na.rm =  T))

sep_same_linear_grades <- 
  general_reg_separe_race_het("ssi_same_negro_na_2_std", 
                              "grades_score", 
                              hetero_var = "grades_score_dummy")
sep_other_linear_grades <- 
  general_reg_separe_race_het("ssi_other_race_negro_na_2_std",
                              "grades_score", 
                              hetero_var = "grades_score_dummy")
sep_all_linear_grades <- 
  general_reg_separe_race_het("ssi_all_race_2_std",
                              "grades_score",
                              hetero_var = "grades_score_dummy")

list_regs_c_0 <- list(sep_all_linear_grades[[1]][[1]][[3]][[7]], 
                      sep_all_linear_grades[[1]][[4]][[3]][[7]],
                      sep_same_linear_grades[[1]][[1]][[3]][[7]], 
                      sep_same_linear_grades[[1]][[4]][[3]][[7]], 
                    sep_other_linear_grades[[1]][[1]][[3]][[7]],  
                    sep_other_linear_grades[[1]][[4]][[3]][[7]])

list_regs_c_1 <- list(sep_all_linear_grades[[2]][[1]][[3]][[7]], 
                      sep_all_linear_grades[[2]][[4]][[3]][[7]],
                      sep_same_linear_grades[[2]][[1]][[3]][[7]], 
                      sep_same_linear_grades[[2]][[4]][[3]][[7]], 
                    sep_other_linear_grades[[2]][[1]][[3]][[7]],   
                    sep_other_linear_grades[[2]][[4]][[3]][[7]])

list_regs_0 <- list(sep_all_linear_grades[[1]][[1]][[1]][[7]], 
                    sep_all_linear_grades[[1]][[4]][[1]][[7]],
                    sep_same_linear_grades[[1]][[1]][[1]][[7]], 
                    sep_same_linear_grades[[1]][[4]][[1]][[7]], 
                    sep_other_linear_grades[[1]][[1]][[1]][[7]],   
                    sep_other_linear_grades[[1]][[4]][[1]][[7]])

list_regs_1 <- list(sep_all_linear_grades[[2]][[1]][[1]][[7]], 
                    sep_all_linear_grades[[2]][[4]][[1]][[7]],
                    sep_same_linear_grades[[2]][[1]][[1]][[7]], 
                    sep_same_linear_grades[[2]][[4]][[1]][[7]], 
                    sep_other_linear_grades[[2]][[1]][[1]][[7]],  
                    sep_other_linear_grades[[2]][[4]][[1]][[7]])

n_obs_0 <- sapply(list_regs_0, get_obs)
n_obs_1 <- sapply(list_regs_1, get_obs)

stargazer(list_regs_c_0,
          column.separate = c(2,2,2),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("all", "rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White", 
                             "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_0)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_grades_0.tex")

stargazer(list_regs_c_1,
          column.separate = c(2,2,2),
          dep.var.labels.include = F,
          dep.var.caption = "Social status",
          column.labels = c("All races", "Same race", "Other race"),
          omit = variables_to_omit, 
          covariate.labels  = c("Grades"),
          omit.stat = c("rsq","ser"),
          add.lines = list(c("Racial group", "White", "Nonwhite", "White", 
                             "Nonwhite", "White", "Nonwhite"),
                           c("Observations", n_obs_1)),
          font.size = "footnotesize", type = "latex",
          out = "tables/hetero/sep_het_lin_grades_1.tex")
